#include <stdlib.h> /* rand & malloc */
#include <stdio.h>  /* printf */
#include <assert.h>
#include <sys/time.h>
#include <cuda.h>

// Integer ceiling division
#define DIV_CEIL(a,b) (a+b-1)/b
#define MICROSEC (1E-6)

#define EPSILON 0.05f

int timeval_subtract (struct timeval *result, struct timeval *x, struct timeval *y);

// Device code
__global__ void fof(unsigned long int N, uint1 *hits, float3 *pts)
{
#define KERNEL_FLOP 6
#define KERNEL_MOP 5
	float d;
	float3 pt;
	// From nbody_kernel.cu
	int i = __mul24(blockIdx.x,blockDim.x) + threadIdx.x;;
	if (i < N) {
		pt = pts[i];
		d = sqrt(pt.x*pt.x + pt.y*pt.y + pt.z*pt.z);
		if (d<EPSILON) {
			atomicAdd(&(hits->x),1); //atomicInc
		}
	}
}

// Host code
int main(int argc, char** argv)
{
	unsigned long int N = 0;
	unsigned int threadsPerBlock = 0;
	unsigned int blocksPerGrid = 0;
	
	// Handle input params
	if (argc<3 || argc>4) {
		printf("Usage: %s 2^<bits> <threadsPerBlock> [<blocksPerGrid>]\n", argv[0]);
		exit(1);
	}
	assert(argc==3 || argc==4);
	int i;
	for(i=0; i<argc; i++) {
		assert(argv[i]);
	}
	N = 1 << atoi(argv[1]);
	threadsPerBlock = atoi(argv[2]);
	if (argc==4) {
		blocksPerGrid = atoi(argv[3]);
	} else {
		blocksPerGrid = DIV_CEIL(N,threadsPerBlock);
	}

	// Hardware limits
	assert(threadsPerBlock<=512);
	assert(blocksPerGrid<=65535); //??
	assert(blocksPerGrid*threadsPerBlock>=N); // enough threads

	// Allocate vectors in device memory
	size_t size = N * sizeof(float3);
	assert(0<size);
	
	float3 *h_pts = NULL;
	float3 *d_pts = NULL;
	uint1 *h_hits = NULL;
	uint1 *d_hits = NULL;

	// Alloc and set host and device memory
	h_pts = (float3 *)malloc(size);
	assert(h_pts);
	h_hits= (uint1 *)malloc(sizeof(uint1));
	assert(h_hits);
	cudaMalloc((void **) &d_pts, size);
	assert(d_pts);
	cudaMalloc((void **) &d_hits, sizeof(uint1));
	assert(d_pts);
	srand48(0);
	for (i=0; i<N; i++) {
		h_pts[i] = make_float3((float)drand48(), (float)drand48(), (float)drand48());
	}
	*h_hits = make_uint1(0);
	assert(h_hits->x==0);
	
	// Copy vector from host memory to device memory
	// h_pts is an input vectors stored in host memory
	cudaMemcpy(d_pts, h_pts, size, cudaMemcpyHostToDevice);
	cudaMemcpy(d_hits, h_hits, sizeof(uint1), cudaMemcpyHostToDevice);

	// Start timer
	double secs = 0.0;
	struct timeval start = {0L,0L}, end = {0L,0L}, elapsed = {0L,0L};
	assert(gettimeofday(&start,NULL)==0);

	// Invoke kernel
	fof <<<blocksPerGrid, threadsPerBlock>>> (N, d_hits, d_pts);
	cudaThreadSynchronize();

	// Stop timer
	assert(gettimeofday(&end,NULL)==0);

	// Copy result from device memory to host memory
	// h_C contains the result in host memory
	cudaMemcpy(h_hits, d_hits, sizeof(uint1), cudaMemcpyDeviceToHost);

	// Check for correctness
	//assert(0 <= h_hits->x);
		
	// Print results
	assert(timeval_subtract(&elapsed, &end, &start)==0);
	secs = (double)elapsed.tv_sec + ((double)elapsed.tv_usec*MICROSEC);
	printf("Secs: %f\n",secs);
	printf("Mops: %f\n",(N*sizeof(float)*KERNEL_MOP / secs)*1e-9);
	printf("Flops: %f\n",(N*KERNEL_FLOP / secs)*1e-9);
	printf("Hits: %d\n",h_hits->x);

	// Free device/host memory
	cudaFree(d_pts); d_pts = NULL;
	cudaFree(d_hits); d_hits = NULL;
	free(h_pts); h_pts = NULL;
	free(h_hits); h_hits = NULL;

	return 0;
}


/*
 * http://www.gnu.org/software/libtool/manual/libc/Elapsed-Time.html
 * Subtract the `struct timeval' values X and Y,
 * storing the result in RESULT.
 * Return 1 if the difference is negative, otherwise 0.
 */

int timeval_subtract (struct timeval *result, struct timeval *x, struct timeval *y) {
	/* Perform the carry for the later subtraction by updating y. */
	if (x->tv_usec < y->tv_usec) {
		int nsec = (y->tv_usec - x->tv_usec) / 1000000 + 1;
		y->tv_usec -= 1000000 * nsec;
		y->tv_sec += nsec;
	}
	if (x->tv_usec - y->tv_usec > 1000000) {
		int nsec = (x->tv_usec - y->tv_usec) / 1000000;
		y->tv_usec += 1000000 * nsec;
		y->tv_sec -= nsec;
	}

	/* Compute the time remaining to wait. tv_usec is certainly positive. */
	result->tv_sec = x->tv_sec - y->tv_sec;
	result->tv_usec = x->tv_usec - y->tv_usec;

	/* Return 1 if result is negative. */
	return x->tv_sec < y->tv_sec;
}

